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Abstract 



We have investigated the relaxational dynamics for a protein model at 
various temperatures. Theoretical analysis of this model in conjunction with 
numerical simulations suggests several relaxation regimes, including a single 
exponential, a power law and a logarithmic time dependence. Even though 
a stretched exponential form gives a good fit to the simulation results in the 
crossover regime between a single exponential and a power law decay, we have 
not been able to directly deduce this form from the theoretical analysis. 
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I. INTRODUCTION 



The manner in which proteins relax to their equilibrium or folded configurations is an 
important subject of current research [I]. In particular the folding kinetics, namely the 
asymptotics of the relaxation of a physical quantity such as the total energy, provides useful 
clues to the general folding problem. In this work we shall provide both numerical and theo- 
retical results concerning the asymptotic behavior of relaxation of a model protein following 
a rapid temperature quench from high to low. We emphasize and attempt to answer the very 
relevant question: what is the true (possible) asymptotics of such a relaxation of a protein ? 
Due to the complex nature of this question, we shall use a simplified protein model for our 
numerical simulations which is based on the monomer-monomer interaction matrix of Ref. 
||; and our analytical work is based on the hierarchically constrained dynamics for glassy 
relaxation ||. 

There existed studies of the relaxational dynamics or related equilibrium properties of 
protein-like heteropolymers where part or the whole of the monomer-monomer interaction is 
modeled by a stochastic term |],|]-|9| . This random interaction clearly gives rise to a compli- 
cated free energy landscape with a large number of highly degenerate energy levels. These 
heteopolymers are believed to simulate certain protein behavior: interesting and important 
results for our intuitive understanding of the folding kinetics have been obtained through 
these studies. Using a lattice copolymer model with Metropolis Monte Carlo dynamics, Ref. 
H examined the folding kinetics as a function of the simulation time t for various tempera- 
ture T, starting from an unfolded initial copolymer configuration. Importantly, the folding 
of the copolymer was found to take place in a two-stage process: a rapid collapse followed 
by a slower adjustment toward the ground state, for a broad range of temperatures. A 
kinetically defined glass transition temperature T g and a thermodynamically defined fold- 
ing temperature 7/ Q ^ were naturally found from the copolymer simulation data ||. These 
temperature scales were used to characterize the foldability of the polymers ||. Refs. |||9| 
directly measured the relaxation of the total system energy E for a random heteropolymer 
after a temperature quench where the interaction parameters between the monomers were 
drawn from a Gaussian distribution. The Monte Carlo M and molecular dynamics [[J data 
were well fit to a stretched exponential form as a function of the simulation time t. One 
additional interesting finding of these simulations was that for longer heteropolymers the 
collapse followed a two-stage process |{J characterized by different sketching exponents, while 
for short chains only one stage was found @. 

Along another line of development, fruitful results have been obtained from analyzing 
even simpler discrete models such as the random energy model [10|. The application of this 
model to protein folding was carried out in Ref. [II] where the folding time was analyzed. 
Interestingly, the random energy model gives different relaxation asymptotics depending on 
the details of the transition matrix between different microstates |12| even though the de- 
tailed balance condition has been enforced. Indeed, power-law [12] decay or logarithmic [II 



decay were both obtained. This is to be compared with the possible stretched exponential 



2 



decay |||9| form. 

Previous investigations indicated clearly that a detailed analysis of the energy level distri- 
bution is one requirement for understanding the related relaxation modes. However, this is 
not sufficient to determine the relaxation behavior completely. It is also important to know 
how the random heteropolymer or protein can make transitions between different energy 
levels corresponding to different polymer configurations. If such transitions are possible, 
one says that the levels are connected and the numerical values of the transition rates will 
determine the ease of transition. However, due to the microscopic origin of the connected- 
ness and the transition rules, it becomes difficult to account for these parameters correctly. 
It is therefore a difficulty problem to predict the true asymptotics of a relaxation in these 
systems. 

From the Master Equation approach point of view, there are difficulties for the investi- 
gation of dynamical systems involving random interactions, and in particular to predict the 
precise form of the relaxation process for folding. Let's consider the relaxation to equilib- 
rium in the unbiased random energy model [14]]. Here it was believed that the relaxation 
to equilibrium due to a temperature quench from above to below the freezing temperature 
follows a two stage kinetics |H |. The first stage is characterized by a logarithmic relaxation 
of the extensive part of the energy, while in the second stage, when the system has already 
reached its frozen state, the relaxation follows a power law decay. However it has also been 
suggested that the relaxation dynamics can be described by a stretched exponential time 
dependence [13], or a single power law decay fll2 |. In all these references a master equation 
was used to investigate the kinetics of relaxation in the context of the random energy model. 

Why are there several different predictions for the asymptotics of the relaxational dy- 
namics for the same model ? The answer is related to the fact that different transition rates 
used in the Master Equation can lead to different asymptotics [ F2| . If we define Wij to be a 
transition rate from the i-th energy level to the j-th one, and define P^ q , Pj q to be the equi- 
librium probabilities for the two states, then according to detailed balance Wy-Ff 9 = WjiPj 9 . 
In terms of the energies of the i-th and j-th levels, the detailed balance condition takes the 
form exp — ^ = Hjjexp — Koper and Hilhorst suggested [12|] a general form of the 
transition rates: Wy = H^exp — (1 — q)^ + q^- where < q < 1. Adjusting the value of q 
may lead to different asymptotics from the Master Equation. Indeed, one finds a stretched 



exponential decay |13| by choosing q = 0, while a power law decay []T4],[T2] for q = 1 (choosing 
q = 1 means that all energy levels have the same activation energy). 

While the relaxational dynamics by quenching to low temperatures is clearly complicated, 
one would in the case of high temperatures intuitively expect a single exponential relaxation 
as can be shown by investigating the high temperature limit using a highly connected random 
master equation. At T — ► oo, detailed balance gives = Wji. Assuming that only a 
fraction 7 of states is connected, a possible choice of the transition rate is = 1 with 
probability 7 and Wy = with probability 1 — 7. A study of the related random master 
equation leads to the conclusion that its eigenvalue distribution is sharply peaked at A^7 
with a width of the order of viV where N is the number of states. This means that there 
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is a dominant time scale which is approximately given by j^-, leading to single exponential 
decay at the high temperature region |L5| . 



Since the random energy model and a highly connected random master equation are 
related to the problem of protein folding (both approaches are combined in the work of Ref. 
||1 1|| ) , we arrive at the following picture for the behavior of the relaxational dynamics of 
proteins. At high temperatures, relaxation is described by a single exponential decay. Upon 
quenching from an unfolded state to equilibrium with a final temperature of the order of 
the freezing temperature, one might expect a power law, or a stretched exponential decay, 
while quenching to a temperature below the freezing point may well lead to an extremely 
slow logarithmic decay. 

In the following sections we examine the above physical picture carefully by means of 
both numerical and analytical methods. Our analysis also provides insight into the crossover 
in the relaxation from a single exponential form at high temperatures to a non-Arrhenius 
relaxation near the freezing temperature. Our numerical simulations are based on Monte 
Carlo dynamics || for a protein model and we fit the data to various possible decay modes. 
In the theoretical model we note the analogy between random heteropolymer relaxation 
and spin glass dynamics, this allows us to construct an analytical model for the protein 
relaxation. The key idea is that the constraints to the relaxation of the degree of freedom at 
each energy level along the relaxation pathway must be taken account 0, and the constraint 
at a lower energy level depends on what happens on the higher one. The asymptotics of our 
analytical model together with the simulation data allows a physically reasonable picture to 
emerge for the relaxational dynamics of proteins. 

The paper is organized as follows. Section II gives the results of numerical simulations for 
the distinct asymptotic relaxation behavior of the folding process at various temperatures 
using random sequences with protein-like energetics 17]. Section III contains our theoretical 
analysis where we derive the asymptotics of the relaxation. Finally Section IV is reserved 
for a short summary. The detailed algebra for the analysis is organized into two Appendices. 



II. NUMERICAL SIMULATION 

We simulated the relaxational dynamics of a model protein by focusing on the energy 
relaxation as a function of the simulation time. The simulation was carried out using the 
Monte-Carlo (MC) method on a 3D cubic lattice for a model protein with 27 amino acid 
residues. The contact interaction between the residues is described by the model of Li, Tang 
and Wingreen ||. In this model, an analysis of the correlations between the elements of the 
Miyazawa-Jernigan (MJ) contact energy matrix [|16|], gave the contact energy of the protein 
in the following form PJT7|, 

E = nq . (1) 

Here n and q are L dimensional vectors. L is the number of residues: L = 27 throughout 
all our numerical simulations. The vector n describes the geometric shape of the protein 
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and its element n« is equal to the number of nearest neighbor contacts made by the i-ih 
residue with other residues, as discussed in Ref. [|I7|. The vector q describes the interaction 
strength of the residues: the i-ih element of q corresponds to the strength of the i-th 
monomer on the chain. The values of qi were found from a fit to the MJ matrix in Ref. [0 . 
Since there are twenty different amino acid residues, there will be twenty different values 
of qi as documented in Ref. 0. As the MJ matrix was obtained from experimental data, 
the elements of q can be viewed as empirical parameters to our model. To use Eq. ([[]), 
we call two residues nearest neighbors if they are not connected along the chain but are 
separated by one lattice constant. In our simulation we fixed the value of q^ in the range of 
—2.5 < qi < 0.0, where a more hydrophobic amino acid residue has a more negative value 
of qi. For example, for the 2D model protein depicted in Fig. (|]), the first integer above 
each monomer corresponds to the number of nearest neighbors for that monomer, while the 
value qi of that monomer is written in parentheses. Then from Eq. ([J), the total energy of 
this 2D protein is E = 3qi + 2q 2 + q% + q± + qe + qs + <?9 + qu + qn- We emphasize again that 
n specifies the geometry of the conformation while q specifies its protein sequence. Finally, 
we comment that in writing down Eq.([l]) we have neglected a quadratic term |2]JT7|] which 
gives a small correction to the linear term of Eq. (|]) . However, the relaxational dynamics 
of the model protein is not substantially affected by this quadratic term. 

Our simulation considered nine randomly generated protein sequences, i.e. nine differ- 
ent q vectors. For each sequence we generated a set of relaxation curves corresponding to 
different quenches to different final temperatures and we studied the decay modes by aver- 
aging over the data for all nine sequences. The allowed moves for a monomer in the Monte 
Carlo simulation were chosen to be the end move, the corner move, and the crankshaft move 
||. For a given protein sequence, the initial state of the model protein was chosen to be 
completely unfolded so that the initial temperature Tj = 16 was about four times higher 
than that of the folding temperature of our model. After the equilibration of the protein at 
the high temperature T i; it was quenched to a set of different final temperatures Tf and the 
energy relaxation curves were recorded. Every quench Tj — > Tf consisted of the following 
operations: a short equilibration of the system at high temperature Tj using 50, 000 MC 
steps followed by a quench to the final temperature Tf. The relaxation data was then stored 
over 200, 000 MC steps. After a relaxation curve was recorded the same procedure was 
repeated many times for thermal averages starting from the same initial condition: 200 to 
1500 quenches were averaged depending on the value of Tf so as to obtain a smooth relax- 
ation curve. Finally, we averaged over a set of different initial equilibrium configurations 
each with the same value of Tj. 

Typical relaxation curves for a particular sequence and a set of different Tf are presented 
in Fig. (0). All curves show a fast decay immediately after the temperature quench followed 
by a slow relaxation to the asymptotic regime. As we were interested in the non-Arrhenius 
relaxation modes, we fitted these relaxation curves by a stretched exponential, a power law 
and a logarithmic form respectively. The stretched exponential decay was taken to have the 
form 
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R(t) = a + a ie xp(-—y . (2) 
a 2 



The power law decay was chosen to be 



R(t) = a + (^r . (3) 

Finally, the logarithmic decay mode was fitted by 

R(t) = a — ailnt . (4) 

In equations a , ai, a 2 , 7, and f3 are fitting parameters. In the following we analyze 

these fittings. 

The fitted value of the stretched exponent 7 is shown in Fig. (|3]) as a function of the 
final temperature Tf. This plot is obtained from a particular sequence q. In the fit we 
used the entire range of the relaxation data where t G (2.0, 10.0) (see Fig. (g)). Although 
there are substantial error bars for 7, its temperature dependence appears to have three 
different regions. At high Tf, the data is consistent with 7 w 1. This is reasonable since the 
relaxational dynamics at a high Tf should, as discussed above, be a single exponential decay. 
For an intermediate range of Tf as shown in Fig. fl3|), 7 = 7(T/) reduces steadily with Tf. 
Finally at very low temperatures, Tf G (0.0,2.0), 7 has a value close to 0.3. However, such 
a small value of 7 usually indicates that the decay is actually not in a stretched exponential 
form but has instead a different relaxational behavior such as a power law. 

The result of fitting the relaxation curves to a power law is summarized in Fig. (|j) , where 
the fitted power (3 is shown as a function of Tf. The difficulty in fitting to a power law is 
that a divergence occurs at the quenching time t — 2 (see Fig. (||)). Since we expect the 
power law decay to be a possible asymptotic of the relaxation, we used the simulation data 
for the time period of t G (2.1, 10.0) to avoid the divergence. Also, as the final temperature 
Tf increases, the characteristic relaxation time (parameter a 2 ) decreases substantially and 
it becomes quite difficult to fit the data to a power law accurately. This is the reason why 
the dependence of (3 on Tf shown on Fig. (|^) is only shown below Tf = 3.0. Nevertheless, 
it is clear that this dependence in the range Tf G (2.0, 3.0) is given by roughly a straight 
line, namely (3 ~ Tf — Tf. Extrapolating this line backwards gives Tf ~ 1.5. An important 
feature of Fig. (||) is that as Tf is in the range G (0.0, 2.0), the relaxation become extremely 
slow and the power (3 levels off at a small value j3 ~ 0.3. The fact that both the stretched 
exponent 7 and the power law (3 change their behavior at nearly the same value of 2.0) 
signals that a transition to a frozen state occurs near this temperature. 

We attempted to fit the relaxation data for the low temperature range Tf G (0.0, 2.0) 
with a logarithmic decay. An almost perfect fit was found for all the sequences at T l f° 9 ~ 1.5, 
as shown in Fig. <^). Note that T l f° 9 ~ Tf. We believe that the logarithmic decay mode sets 
in near the freezing temperature as demonstrated by the random energy model ||14|| . We thus 



denote this special temperature scale to be T G : T G ~ T l f° 9 ~ Tj 3 . We shall use the notation 
Tq in the rest of this paper as it is reminiscent of a glass transition temperature. Note 
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that the nature of glass transition temperature for copolymer systems has been thoroughly 
examined in Ref. || and we refer interested readers to that article for details concerning the 
related physics. 

Our numerical prediction of the occurrence of a low temperature logarithmic decay mode 
at temperature ~ Tq agrees with the prediction for the random energy model of Ref. [Ill 



In addition, from the theoretical analysis of the random energy model in Ref. [jEJ, it follows 
that the relaxation near Tq should follow a power law decay with a power /3 ~ (T—Tq). Our 
simulation data presented in Fig. (0!) is consistent with this linear dependence. On the other 
hand, we found that fitting the data near the freezing point Tq with a stretched exponential 
does not produce an exponent 7 lying along the line which would be expected from Ref. Jl3| . 
Finally, from our numerical data it still remains unclear whether the asymptotics above the 
freezing temperature follows a stretched exponential decay, although it is apparent that the 
data can be reasonably well fitted by such a stretched exponential. 



III. THEORETICAL ANALYSIS 

To understand the simulation data presented in the last section, in the following we make 
a theoretical analysis of the relaxation based on the model of hierarchically constrained 
dynamics for glassy system relaxation 0. 

The conventional Arrhenius relaxation is described by a single exponential decay with a 
time scale r 

R{t) = Roexp{--) (5) 

T 

where R(t) is the physical quantity of interest. We obtain a more complicated relaxation 
behavior by assuming a distribution of time scales. Suppose that P(r) is a continuous 
distribution of possible time scales in our system, then, using the notion of parallel relaxation 
we write for the total relaxation process || 

/Tffi ax 
P(r)exp( )dr . (6) 
min T 

This approach is especially attractive in the case where dynamics is modeled by the Master 
Equation. As the Master Equation is a first order differential equation, the general solution 
will be a sum of exponential decays with time scales inversely proportion to the eigenvalues 
A of the transition matrix. A distribution of time scales P(r) thus corresponds to the 
distribution of inverse eigenvalues P(j)- 

As pointed out by Palmer et. al. [|], though a description of relaxation by (|6]) is certainly 
appealing, the nature of P(r) is not a priori clear. In the dynamics of such complex systems 
as random polymers or proteins, it is intuitively plausible that P(r) should take into account 
the dynamic constraints: e.g. monomer A cannot move until monomer B moves out of the 
way. P(t) may also depend on such factors as ergodicity breaking in a frozen state of the 
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system. The theory should also involve a hierarchy of degrees of freedom from fast to slow. 
For the protein problem, the fastest degrees of freedom could involve single-atom motion, 
while slower degrees of freedom could involve a diffusion of domains of secondary structure. 

The original theory of Ref. || considered a discrete set of levels n = 1, 2, 3..., iV with the 
degrees of freedom in level n represented by N n pseudospins 5*. Each spin in the (n — l)th 
level is only free to change its state if a condition on some spins at level n is satisfied, thus 
providing a hierarchy of constraints. This is the key idea of the model as shown schematically 
in Fig. (0). We take the constraint to be that jl n spins in level n (jl n < N n ) attain one 
particular state out of their 2^ n possible states. Then the average relaxation times r n _i will 
be related to r n as r„_i = 2^ n r n . Hence, 

N 

r n = r N exp(J2 A*i) (7) 

i=n 

where ^ = jliln2, and tn sets the time scale of the relaxation. The average relaxation can 
now be written as follows 

N t 

R(t) = RoY, N n™P( ) (8) 

71=1 Tn 

N 

T n = T N exp(%2 fXi) . (9) 

i=n 

To apply this theory to protein folding dynamics, we notice that the summation is 
taken over (N — n) higher energy levels in Eq. (|7|) since the states corresponding to the 
higher energy levels (the unfolded states) relax very fast because of minimal dynamical 
constraints. Clearly, as the folding proceeds and the level index n decreases, there will be an 
increased number of dynamical constraints on folding as the protein becomes geometrically 
more compact. Hence the relaxation times for the more folded states should increase in 
comparison to those of the less folded states. The meaning of the "pseudospins" of Ref. 
|| in our protein problem can be elucidated by noticing the similarity of a spin flip and a 
local monomer move. For instance, a crankshaft move, a corner move or an end move in our 
Monte Carlo simulation could be compared with a pseudospin flip, while a configuration of 
spins could represent a protein conformation. Two energetically close protein conformations 
may differ from each other by a local move of a monomer, similar to the flip of a spin. It 
is much easier to make a local move (a spin flip) in an unfolded state than in a folded state 
because of the energetics consideration and the importance of volume exclusion interactions 
in compact protein conformations. Clearly, the analogy of pseudospin flips to local monomer 
moves can be extended to more realistic off-lattice protein models. 

From the above discussion, we now make some assumptions about the form of /ij and 
N n for the protein problem. First, from the form of Eq. (H) one can immediately deduce 
that fa should depend upon temperature. When T — > oo, the relaxation should be a single 
exponential, from Eq. (||) we conclude r n (T — > oo) = const and it should not depend 
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on n. The simplest form for «j satisfying these conditions is «j = -jr. As far as N n is 
concerned, it should be proportional to the total number of states on each level. When 
T — > oo, all the pseudospins are excited so that N n in this limit equals to the number of 
pseudospins on the nth level. As temperature decreases and the system finally goes into 
the frozen state, the number of pseudospins on each level should decrease to the order of 
unity signifying ergodicity breaking in the relaxational dynamics. For our protein problem, 
in the low temperature limit the decrease in the number of pseudospins corresponds to the 
restriction of the dynamical space to a set of specific pathways for protein relaxation. In our 
model we choose N n to depend on temperature in the form 

N n = A" (1 ~ Ms) (10) 

where u g = and corresponds to the "level freezing" temperature of Ref. Sum- 
marizing all these considerations, our complete model for protein relaxational dynamics is 
given by 

R(t) = i?o E A n(1 ~^W( ) (11) 

71=1 T " 



N rpf 

Tn = T N eX P (j:^) . (12) 

i=n 

It is convenient to write our model in an integral form, so that 

r N t 
R{t) = Ro N{n)exp( )dn . (13) 

Jl T n 

Here the function N(n) is the continuous version of fllQf) . 

To make further analysis tractable and considering the fact that iV is large, we make a 
reasonable assumption that and T£ to be roughly constants for all energy levels, so that 
fig = ?y and fif = ^- . The above expressions can be written in terms of the distribution of 
time scales r: 

R(t) = R [ Tma \N(nA\ T=Tn ex P (--)dT (14) 

J T min (XT T 

where r m j n = and r max = r^exp (X)i=i ^i)- Comparing the above equation with @, we 
notice that the distribution of time scales is just 

P(r) = N(n)^\ T=Tn , (15) 

where the right hand side is understood as replacing all n dependence by r through the 
transformation r = r n using Eq. (|T2|). It is then straightforward to obtain, 



lax ~t 

R{t) = R I P(r)exp(--)dT (16) 
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with an unnormalized distribution of time scales 

Tat a-Mg^A 

P(t)~(-) "t ■ (17) 

T 

Denoting 

= ^^InX = T^-ln\ (18) 
we obtain a main result of this analysis, 

( — f +1 exp(--)dr . (19) 

- min T T 

One of the important outcomes of our model is that the distribution of time scales is a power 
law. We now discuss several possible cases for different temperatures. 



A. The case of T < T 9 

For low temperatures T < T 9 , from Eq. (|T^) /3 < 0. We define rj = —[3 so that rj > 0. 
Then the normalized time scale distribution ([T7|) becomes 

l-r, 



P(r) = - 71 



( Tmax \ 
Train J 



(20) 



The relaxation, obtained from Eq. (19), is 



Tmax 7~min 



R{t) = v ' v / r" X- l - v exp(-x)d X (21) 



where T min = and r max = T^exp{N^r). R(t) is normalized so that R(t)\ t =o = 1-0. Notice 
that we cannot simply extend the integration from (r min , r max ) to (0, oo) because the integral 
would then diverge. Still, it is not difficult to obtain an asymptotic form of this expression 
as i] — > and t ^> T m j„. From the detailed analysis presented in Appendix I, the final 
asymptotic of the relaxation for this temperature range is given by, 

R(t) x ci - c 2 lnt (22) 

where c\ and c 2 are constants. 

Fitting the low temperature simulation data of the last section to this logarithmic decay, 
as shown in Fig. (|5|), a very good fit is obtained in the interval 1.5 > T > 1.25. Since T 9 is the 
level freezing temperature, and when all the individual energy levels freeze, the whole system 
freezes. Thus we must have T 9 ~ Tq which is the freezing transition temperature discussed 
in section II. This allows us to conclude that the hierarchically constrained dynamics model 
presented in this section has a logarithmic decay mode for low temperatures T < Tq as the 
true asymptotic. 
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B. The case of T > T 9 



When at higher temperatures T > T 9 , from Eq. fll8|) we have > 0. In this Case, cL 
normalized time scale distribution, QI7|), will be 



P(r) 









(t ■ \P~ 


1 - 


1 min \ 




\ Tmax / 



1+/3 



(23) 



The relaxation is thus 
J2(f) 







T rmn 'max 



' mtn 

t 



1 exp(-x)dx ■ 



(24) 



As r ma;r ;§> t ^> r m j n , the integral in can be extended from (r min , r max ) to (0, oo) without 
introducing too much error (the integral is also convergent in the whole interval). Hence, 



R(t) x const(/3)( 



t 



(25) 



We thus obtain a power law decay for the relaxation when T > T 9 . Moreover, since the 
exponent (3 = T ~J 3 InX, we can expect a linear dependence of (3 as a function of temperature 
near the freezing transition point. Our simulation data of Fig. (Q) agrees with this result 
and gives T 9 = 1.5, ££■ = 1.2 . 



C. The case of T > T 9 

We note that formally, the case of (3 > should be divided into two regimes: 1 > > 
and /? > 1. The latter has a special point = 1 because of the integrand of Eq. ( p4j ) which 
can be written as follows 

e(j3,t)= f~ x^expi-xW ■ (26) 
Jo 

Taking into account the continuous behavior of Q(0,t) for > 1, let us calculate the 
relaxation 

R(t) = m^m. ^, t) (27) 

as = p where p is any positive integer. In this case 0(0, t) can be calculated analytically 
and the details are given in Appendix II. The relaxation function is found to be 

RW ^ + jM) + (P IM+ 2) + (p + + 3) + -» (28) 

Hence, as temperature increases, non-Arrhenius decay modes would be observed as asymp- 
totics at t ^> (p + l)r min , while an initial stage of relaxation t < (p + l)r min will be an 
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Arrhenius-like single exponential decay. As temperature increases, the overall relaxation 
time to equilibrium becomes shorter and, finally, it becomes impossible to distinguish non- 
Arrhenius asymptotics at the tail of the relaxation function. Notice that the onset of the 
single exponential decay occurs at temperatures corresponding to (3 <^ 1 or T <^ T 9 + f^-. 
From the simulation data of Fig. (|4]), we can calculate that this occurs as T > 2.7 which 
corresponds approximately to the midpoint of the crossover regime from 7 = 1.0 to 7 ~ 0.3 
on Fig. ([|). This result indicates that although one can fit the simulation data as a stretched 
exponential, it is quite possible that this fit does not in fact give the correct asymptotic na- 
ture of the crossover from a single exponential to a power law decay as temperature changes. 

Finally, we comment that we were not able to obtain a stretched exponential asymptotic 
decay from our analytical theory. In general, a transition from single exponential decay at 
high temperatures to power law decay at lower temperatures could proceed via a set of more 
intricate functions reminiscent of stretched exponentials. However, whether these functions 
have stretched exponentials as their true asymptotics can not as yet be resolved from the 
theory presented in this section. 

IV. SUMMARY 

In this work we addressed the problem of the relaxational dynamics of a model protein 
for temperature quenches from an unfolded state at high temperature to a range of low 
temperatures. The results of our Monte Carlo simulations and our theoretical analysis of 
a hierarchical model proposed for the protein relaxation dynamics clearly indicate three 
different decay modes: first a logarithmic decay at T < T 9 , next an asymptotic power law 
decay at T > T 9 with the power (3 ~ (T — T 9 ) and finally, as the temperature increases to 
T > T 9 + p-, a single exponential decay. Here the temperature scale T 9 describes a freezing 
transition where the decay mode changes. Our numerical and analytical results suggest that 
a crossover from a power law relaxation to a single exponential decay occurs via an intricate 
interplay between the time intervals for which these decay modes are valid, although this 
crossover can be reasonably fitted numerically to a stretched exponential form. 
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APPENDIX I 



In the appendix we derive the logarithmic decay form for low temperatures T <T g . As 
stated in section III, the relaxation function at (3 < 0, n = —f3 is 



R(t) 



Tjt 71 



Tmax 



mm r max 

We can rewrite this expression in a form 

t 



__ I r ^ 1 V exp (_ x ^ dx _ 



R(t) 



max mm Tmax 



00 f— IV 

,V, J- i = Q 1- 



After integration of the \ variable, we have 

V 



R(t) = 1 + — 



g (-1)7' / I 



Tmax T^in |_ i=1 1]) \T~ mi n Tmax, 

As n — > and noticing T max ^> T min , R(t) becomes 



R(t) = 1 + 



ln- 



E 



-iy ( t 



l\l \T, 



+ 0( V ) 



To sum up the series, we differentiate R(t) to obtain 



R(t)> 



t 



i-l 



.i=l 



\ ' mm / 

We can now perform the summation to obtain 

R(t)' = — — ^ — r— e~ T ™» - 1 

l n (lmax\ t 

As £ ^> r m j n , -R(t)' becomes even simpler 



+ 0(77) 



/n (I***-) t 



+ 0{n) 



Final integration yields 



R(t) = const — 



In (^^) 



ln(t) + 0{n) 



(29) 



(30) 



(31) 



(32) 



(33) 



(34) 



(35) 



(36) 



which is the logarithmic decay. 
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APPENDIX II 

For (3>0, 

t 

Q(/3,t)= ( Tmm x^exp-xdx (37) 
J o 

The relaxation function is then given by 

fl (*) =3g pCTTpCT ^eCg.*) ■ (38) 
Choosing (3 — p, where p is a positive integer, and assuming that r max ^> T min , we obtain 



®(M = Tm * n X p ~ l exp{-x)d X (39) 

*/ 

and 

i2(t)=p(^)*e(p,t) • (40) 
0(/3, t) is a standard integral and can be evaluated to give 

e(/?,*) = 3£>(p-l)!(l-ezp( )(1 + + ^f- + ... + f mm _ )) . (41) 



Noticing that 



1 + ^T + ^^ + - + ¥^T)T" exp( w ) "^^^ (42) 



after substitution to ([0]) and QiOD we find 



t \ ( t \2 ( _L_ \3 



m = exp{ -^- ){1 + (iTr) + ¥Tm^) + (P+ i)(^2) (P+ 3) + ■■•> < 43) 

which is the desired result. 
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FIGURES 

FIG. 1. A 2D lattice model of a protein. The figure over each monomer indicates the number of 
its nearest neighbors. The parameter in parentheses indicates the strength of amino-acid situated 
on a particular site. 

FIG. 2. A typical set of relaxation curves for a temperature quench Tj — > Tf, for a given 
sequence and different final temperatures. The unit of time is 25, 000 Monte Carlo steps. 

FIG. 3. Fit to a stretched exponential decay. All nine 7 — Tf curves for different sequences 
were averaged to obtain the averaged behavior of the stretched exponential relaxation as a function 
of the final temperature of the quench. An individual 7 — Tf curve for a particular sequence was 
obtained by fitting the relaxation curves to a stretched exponential form. 

FIG. 4. Fit to a power law decay. All nine (3 — Tf curves for different sequences were averaged 
to obtain the average behavior of the power law relaxation as a function of the final temperature 
of the quench. An individual (3 — Tf curve for a particular sequence was obtained by fitting the 
relaxation curves to a power law decay form. 

FIG. 5. Fit to a logarithmic decay mode for a specific sequence. Similar fits for other sequences 
show that the best fit is found for Tf ~ 1.5. 

FIG. 6. Model of hierarchically constrained dynamics for relaxation. The probability of the 
pseudospin flip on the (n — l)th level is coupled to the probability of occurrence of a specific 
configuration of pseudospins on level n. 
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